Introducción al manejo de datos geográficos¶
Clase 2¶
import pandas as pd
import numpy as np
import geopandas as gpd
from shapely.geometry import Point
import pyproj
fachadas_sndg = pd.read_csv('../data/fachadas_sndg.csv')
fachadas_sndg.head(2)
fachadas_usig = pd.read_csv('../data/fachadas_usig.csv')
fachadas_usig.head(2)
def construye_coords(x,y):
try:
geometria = Point(float(x), float(y))
return geometria
except:
return 'sin coordenadas'
def construye_gdf(df,proj):
geom = df.apply(lambda x: construye_coords(x.lon, x.lat), axis=1)
df['geometry'] = geom
localizables = df.loc[df['geometry']!='sin coordenadas']
if type(proj) == str:
crs = pyproj.CRS.from_user_input(proj)
else:
crs = pyproj.CRS(proj)
gdf = gpd.GeoDataFrame(localizables, crs=crs, geometry=localizables.geometry)
print('Devolviendo geodataframe con {} casos perdidos'.format(len(df)-len(gdf)))
return gdf
wkt = """PROJCS["GKBA",
GEOGCS["International 1909 (Hayford)",
DATUM["CAI",
SPHEROID["intl",6378388,297]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",-34.6297166],
PARAMETER["central_meridian",-58.4627],
PARAMETER["scale_factor",0.999998],
PARAMETER["false_easting",100000],
PARAMETER["false_northing",100000],
UNIT["Meter",1]]"""
gdf_usig = construye_gdf(fachadas_usig,wkt)
gdf_sndg = construye_gdf(fachadas_sndg,wkt)
Ahora que finalmente tenemos nuestros geodataframes de puntos listos veamos si las dos herramientas nos devuelven los mismos resultados. Para eso, repasemos algunos conceptos básicos de visualzación. Cuando trabajamos con algunas librerías como matplotlib, es necesario tanto la figura donde se van a plotear nuestros gráficos como el lugar que va a ocupar. En líneas generales, lo que nos debe quedar en claro es que la instancia que nosotros creemos para la figura será siempre el método figure del módulo pyplot de matplotlib. Esto nos permitirá definir atributos como el tamaño. Para la posición de la figura, se utiliza el método add_subplot (aplicable siempre a un objeto de tipo pyplot. Así, se define la forma en la que se disponen las figuras. En nuestro caso, si queremos que el primer mapa se disponga al lado del segundo o, si por ejemplo, queremos superponerlos.
import matplotlib.pyplot as plt
%matplotlib inline
# creamos la figura con un tamaño específico
fig = plt.figure(figsize=(15,5))
# disponemos dos ejes en una columna. El tercero que pasamos como parametro es el número de eje
ax1 = fig.add_subplot(1,2,1)
# acá vemos el eje #2
ax2 = fig.add_subplot(1,2,2)
fig = plt.figure(figsize=(15,5))
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)
gdf_usig.plot(ax=ax1, color='red')
gdf_sndg.plot(ax=ax2, color='blue')
ax1.set_axis_off()
ax1.set_title('Resultados Normalizador USIG')
ax2.set_axis_off()
ax2.set_title('Resultados Normalizador SNDG');
# Si quisieramos verlos superpuestos
fig = plt.figure(figsize=(15,5))
ax = fig.add_subplot(1,1,1)
gdf_usig.plot(ax=ax, color='red')
gdf_sndg.plot(ax=ax, color='blue')
ax.set_axis_off()
ax.set_title('Resultados USIG/SNDG');
Como se puede apreciar, los resultados son bastante similares. Con lo cual, a nivel performance ambos normalizadores trabajaron de manera muy parecida. Pero dónde es que caen esos puntos realmente? En los plots anteriores se puede ver con claridad la forma de la ciudad. De todos modos, es muy común cuando trabajamos con este tipo de situaciones en las que ploteamos algo cuyo resultado inicial nos resulta incierto, querer ver nuestros datos con una especificidad mayor. Hacer zoom, ver en qué polígono cae un punto, el nombre de una calle, etc. Para esto, existen varias librerías que permiten, de una manera muy ágil, plotear figuras sobre un layer dinámico (como si fuese un mapa web). A continuación daremos un pantallazo de algunas de ellas...
1. Leaflet
# importamos la librería
import mplleaflet
# Instanciamos el plot
ax = gdf_usig.head(1000).plot(markersize = 50, color = "red")
# Lo visualiamos con el método display
mplleaflet.display(fig=ax.figure)
Si prestamos atención al warning, se nos sugiere utilizar display de Ipython. Para ello, vamos a tener que, primero, exportar nuestro mapa y traer la ruta a ese ħtml que exportamos. Con leaflet también podemos hacer esto.
# Ploteamos la data,
gdf_usig.head(1000).plot(markersize = 50, color = "red");
# la exportamos
mplleaflet.show()
# y la llamamos nuevamente. IFrame nos permite regular el alto y el ancho del output
path = '_map.html'
from IPython.display import IFrame
IFrame(path, width=1000, height=400)
2. Contextily
Otra es Contextily, una librería muy útil para plotear puntos, polígonos o líneas sobre un layer estático que nos de un primer pantallazo del mapa que buscamos construir.
# importamos la librería
import contextily as ctx
puntos = gdf_usig.plot(figsize=(10, 10), alpha=0.5, color='red', edgecolor='k')
ctx.add_basemap(puntos,zoom=12, crs=4326)
puntos.set_axis_off()
Recuerdan que en la clase anterior, vimos que en la nueva release de Geopandas el sistema de proyección de coordenadas (o CRS) había incorporado una mejora en su manera de definirlo? Y recordemos también que, el motivo de este cambio se apoya en la idea de sustitur un proj4string (o simplemente un string con información sobre el sistema de proyección utilizado) por un nuevo objeto enriquecido que nos permite acceder a distintos tipos de atributos y métodos. Por ejemplo:
# esta es la información de nuestro crs
gdf_usig.crs
# con la que podemos acceder, por ejemplo al nombre
gdf_usig.crs.name
# o al datum
gdf_usig.crs.datum
Sin embargo, esta nueva forma de definir el objeto CRS puede tener un impacto importante en el uso de distintas librerías. Y seguramente llevará un tiempo hasta que esto termine siendo totalmente compatible o todas se adapten al uso de este nuevo objecto. El principal causal de esto que existen muchas fuentes y formas de definir un CRS, algunas de las cuales pueden tener una descripción que no se ajusta completamente a los nuevos estándares de PROJ> 6 (cadenas de proj4, formatos WKT más antiguos, etc.). Como se menciona en la documentación oficial de geopandas sobre el uso de proyecciones, en estos casos, el objeto pyproj.CRS que se obteiene, podría contener algunas incosistencias. Por ejemplo, un código EPSG que no se condice con la localización esperada.En el blog de Joris Van der Boschee se brindan algunas explicaciones adicionales bastante interesantes.
Pero por qué mencionamos esto ahora? Si prestan atención, en la primera instancia del mapa que construimos con Contextily el crs que definimos fue el clásico 4326. ¿Por qué hicimos esto si nosotros construimos el gdf a partir de un wkt? Veamos qué hubiese pasado si definíamos el CRS de otra manera.
# Por ejemplo, si hubiésemos querido asignar el crs de otro dataframe cuyo epsg nos es familiar
barrios = gpd.read_file('../carto/barrios_badata.shp')
barrios.crs.datum
# Primero vemos que nos devuelve un error relativo al crs que definimos desde el wkt
puntos = gdf_usig.to_crs(barrios.crs).plot(figsize=(10, 10), alpha=0.5, color='red', edgecolor='k')
ctx.add_basemap(puntos,zoom=12, crs=barrios.crs)
# Creamos nuestro gdf pero ahora no intentando reproyectar sino asigando el CRS deseado directamente
gdf_usig_crs_barrios = construye_gdf(fachadas_usig,barrios.crs)
# Y volvemos a intentar la operación
puntos = gdf_usig_crs_barrios.plot(figsize=(10, 10), alpha=0.5, color='red', edgecolor='k')
ctx.add_basemap(puntos,zoom=12, crs=gdf_usig_crs_barrios.crs)
Contextily sigue sin reconocer nuestro sistema de coordenadas. Intentemos un camino diferente...
# Veamos en la descripción de nuestro objeto CRS cual es el código EPSG asignado
gdf_usig_crs_barrios.crs.datum
# Y intentemos nuevamente...
puntos = gdf_usig_crs_barrios.plot(figsize=(10, 10), alpha=0.5, color='red', edgecolor='k')
ctx.add_basemap(puntos,zoom=12, crs=6221)
Pueden intentarlo tanto con este código como con el que vimos inicialmente en la descripción de nuestro primer CRS (9001). Ambos devolverán un mensaje de error de CRS no soportado. Esto, porque Contextily no está reconociendo la definición de ninguno de los dos (el que heredamos de nuestro shape de barrios o el especificado en el wkt).
Y qué hubiese pasado si apelábamos al EPSG oficial que corresponde a la faja donde cae la Ciudad de Buenos Aires? Pruébenlo y cuéntennos cómo les va!
De manera diferente, esto no nos hubiese pasado, por ejemplo, con leaflet. Como dijimos antes, el impacto que tiene este nuevo objeto no es siempre el mismo. Veamoslo:
# Instanciamos el plot con nuestro nuevo gdf
ax = gdf_usig_crs_barrios.head(1000).plot(markersize = 50, color = "red")
# Y lo visualiamos con el método display
mplleaflet.display(fig=ax.figure)
# Creamos el gdf con el EPSG que vinos inicialmente en datum
gdf_usig_9001 = construye_gdf(fachadas_usig,9001)
# Y repetimos la operatioria
ax = gdf_usig_9001.head(1000).plot(markersize = 50, color = "red")
mplleaflet.display(fig=ax.figure)
# Este sería el resultado con Contextily si no se especificara el CRS
puntos = gdf_usig_crs_barrios.plot(figsize=(10, 10), alpha=0.5, color='red', edgecolor='k')
ctx.add_basemap(puntos,zoom=12)
En este tipo de circunstancias, algo que es muy recomendable es trabajar directamente en 4326 (para asegurarnos que no tendremos ningún problema de incompatibilidad o falta de soporte). O bien, este nuevo objecto pyproj.CRS cuenta con un método to_epsg() para devolver un id equivalente. Este es bastante útil para identificar si nuestra proyección (en nuestro caso definida desde un wkt) está siendo identificada correctamente...
# En nuestro caso, vemos que el epsg es un None type, algo que ya nos debería parecer sospechoso
print(gdf_usig.crs.to_epsg())
# También podemos acceder a la descripción original si hacemos:
gdf_usig.crs.source_crs
# Y veamos también que ahora trabajamos con un nuevo parámetro: min_confidence
gdf_usig.crs.source_crs.to_epsg(min_confidence=20)
Esto se encuentra bastante bien explicado en la documentación de pyproj. Si el método to_epsg no nos devuelve ningún código EPSG que coincida con el wkt o el pyproj4 string que usamos inicialmente (porque no lo encuentra o porque no existe), podemos apelar al parámetro min_confidence. Este, nos permitiría obtener el código EPSG más cercano, alternando los umbrales que estamos dispuestos a tolerar.
# Vemos que ahora Contextily interpeta el crs...
puntos = gdf_usig.plot(figsize=(10, 10), alpha=0.5, color='red', edgecolor='k')
ctx.add_basemap(puntos,zoom=12, crs=gdf_usig.crs.source_crs.to_epsg(min_confidence=20))
En este caso, si bien el resultado es bastante similar al que obtuvimos con un EPSG igual a 4326, este método nos puede resultar de utilidad para situaciones en las que necesitemos trabajar con una proyección diferente asegurándonos que el mismo sea interpretable.
3. Folium
Folium es una librería que combina distintos tipos de funcionalidades. Desde controles de zoom y la posibilidad de alternar layers hasta clases que devuelven mapas de calor y coropletas.
# Si no lo incluiste en los requirements, lo podemos instalar directo desde el notebook
#!pip install folium
# Importamos la librería
import folium
# Creamos un nuevo gdf en 4326 para los resultados de ambos nomralizadores
sndg_4326 = construye_gdf(fachadas_sndg,4326)
usig_4326 = construye_gdf(fachadas_usig,4326)
# Convertimos nuestro gdf de puntos en geojson
sndg_gjson = folium.features.GeoJson(sndg_4326.head(1000), name="SNDG")
# Acá podemos ver la estructura
#points_gjson.data.get('features')
# Y ploteamos el mapa...
mapa = folium.Map(location=[-34.6157437, -58.4333812],
zoom_start=11,
control_scale=True)
sndg_gjson.add_to(mapa);
mapa
def mapa_folium(gdf, polygons, barrio, exporta):
'''
Construye un mapa de puntos y/o polígonos sobre
un layer dinámico de leaflet.
...
Argumentos:
gdf(GeoDataFrame): Geodataframe de puntos.
barrio (str): nombre del barrio.
polygons (GeoDataFrame): GeoDataFrame de polígonos.
exporta (bool): True para exportar.
Devuelve:
folium.Map : mapa de puntos/polígonos
'''
# Creamos el layer de leaflet. Lo centramos en la Ciudad de Buenos Aires
centroide = gdf.geometry.centroid
coordenadas = [centroide.y.mean(), centroide.x.mean()]
layer = folium.Map(location=coordenadas,
zoom_start=11,
height = 500,
control_scale=True)
# Si pasamos polígonos de base, también los cargamos al mapa
if polygons is not None:
estilo = lambda x: {'fillColor': 'red' if
x['properties']['BARRIO']== barrio else
'grey'}
pol_hover= folium.features.GeoJsonTooltip(
fields=['BARRIO','COMUNA'],
aliases=['Barrio:','Comuna:'])
folium.GeoJson(polygons.to_crs(gdf.crs),
name = 'Barrios de la ciudad',
style_function = estilo,
tooltip = pol_hover
).add_to(layer)
# Agregamos nuestros puntos con el método CircleMarker
for i, row in gdf.iterrows():
folium.CircleMarker((row.lat,row.lon),
radius=3, weight=2,
color='blue',
fill_color='blue',
fill_opacity=.5,
tooltip= 'Certificado:'+row.vencimiento_certificado).add_to(layer)
# Exportar el resultado
if exporta:
layer.save('mapa.html')
folium.LayerControl(autoZIndex=False, collapsed=False).add_to(layer)
return layer
# ploteamos un mapa sólo de puntos
mapa_folium(sndg_4326.head(50), None, None, False)
# ahora lo combinamos con un gdf de polígonos
mapa_folium(sndg_4326.head(1000), barrios, 'VILLA LUGANO', False)
fig = plt.figure(figsize=(15,5))
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)
barrios.plot(ax=ax1, color='grey')
usig_4326.to_crs(barrios.crs).plot(ax=ax1, color='red', markersize=1)
barrios.plot(ax=ax2, color='grey')
sndg_4326.to_crs(barrios.crs).plot(ax=ax2, color='blue', markersize=1)
ax1.set_axis_off()
ax2.set_axis_off()
ax1.set_title('Resultados Normalizador USIG')
ax2.set_title('Resultados Normalizador SNDG');
fig = plt.figure(figsize=(15,5))
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)
barrios.plot(ax=ax1, color='grey')
usig_4326.to_crs(barrios.crs).plot(ax=ax1, color='red', markersize=10, edgecolor='white')
barrios.plot(ax=ax2, color='grey')
sndg_4326.to_crs(barrios.crs).plot(ax=ax2, color='blue', markersize=10, edgecolor='white')
ax1.set_axis_off()
ax2.set_axis_off()
ax1.set_title('Resultados Normalizador USIG')
ax2.set_title('Resultados Normalizador SNDG');
from folium.plugins import HeatMap
def mapa_de_calor(gdf, radio, inicio, fin):
'''
Construye un mapa de calor sobre
un layer dinámico de leaflet.
...
Argumentos:
gdf(GeoDataFrame): Geodataframe de puntos.
radio (int): radio de cada punto.
inicio (int): límite inferior.
fin (int): límite superior.
Devuelve:
folium.Map : mapa de calor
'''
# Ubicamos el centro del mapa
centroide = gdf.geometry.centroid
coordenadas = [centroide.y.mean(), centroide.x.mean()]
# Conseguimos la latitud y la longitud de cada punto con una función anónima
x = gdf["geometry"].apply(lambda punto: punto.x)
y = gdf["geometry"].apply(lambda punto: punto.y)
# Convertimos esas coordenadas en una lista de tuplas
puntos = list(zip(y, x))
# Creamos el layer de leaflet
layer = folium.Map(location=coordenadas,
zoom_start=11, control_scale=True)
# Agregamos el mapa de calor al layer que habíamos instanciado
HeatMap(puntos[inicio:fin],
radius=radio).add_to(layer)
return layer
# también se puede jugar con otros parametros como max_zoom o el tile para elegir otro base map
mapa_de_calor(usig_4326, 15, 1000, 2000)
Para concluir con esta práctica, veamos qué se puede hacer con los datasets que hemos venido construyendo. Un aspecto muy importante en este tipo de situaciones en las que la construcción de nuestros set de datos se vuelve un proceso lento y demandante, es no olvidar o dejar de lado la pregunta que nos va a ayudar a responder. Y no por la pregunta en sí, ya que puede ser una o más. Sino porque, la comunicación de nuestros hallazgos se vuelve mucho más clara cuando respondemos algo en concreto. Por ejemplo, dónde se aglomera nuestra variable de interés? existe algún patrón de localización específico?
Con nuestro mapa de calor, ya pudimos identificar algunas posibles zonas. Sin embargo, vamos a aprovechar esta pregunta para introducir una nueva geometría, un tipo de polígono que no vimos hasta ahora. Este no es nada más y nada menos que el radio censal. Una división o unidad geográfica que utiliza el Instituto Nacional de Estadísticas y Censos (INDEC) para sus relevamientos censales.
Pero no vamos a ahondar mucho en este tema. Simplemente, mencionaremos que una de sus principales ventajas, es la posibilidad de contar con información sociodemográfica a un alto nivel de desagregación territorial (algo que se vuelve importante en términos de acceso a información útil para describir con detalle áreas más reducidas - o con mayor zoom). Tema que retomaremos en otro notebook.
A continuación, veremos un ejemplo en el que utilizaremos los radios solamente por el tamaño de los polígonos (también existen otras librerías muy útiles para lograr este tipo de efectos en nuestros outputs). Haremos un join espacial para ver la cantidad de puntos (nuestras fachadas registradas en las partidas de conservación patrimonial) que caen dentro de nuestros polígonos censales y evaluaremos la concentración de actividad (personas que hayan registrado una fachada) en base a la clasificación de este valor dentro de los radios.
Como vimos en la clase anterior, Geopandas cuenta con un módulo bastante útil que nos permite detectar la relación espacial entre dos geometrías. Este se denomina spatial join y nos permitira conocer cuántos puntos caen dentro del área de cada polígono.
# importamos el módulo desde gepandas
from geopandas import sjoin
# cargamos el shapefile de radios censales
radios = gpd.read_file('../carto/informacion_censal_por_radio_2010.shp')
# mergeamos ambos gdf
fachadas_r = gpd.sjoin(sndg_4326, radios.to_crs(sndg_4326.crs))
# visualizamos el reultado...
fachadas_r.plot();
# Agrupamos la cantidad de fachadas con certificado de conservación por radio censal
fachadas_por_radio = fachadas_r.groupby(['RADIO_I'])[['partida']].count().reset_index()
# Creamos un diccioanrio donde la key es el id del radio y el value la cantidad de partidas.
d = dict(zip(fachadas_por_radio['RADIO_I'], fachadas_por_radio['partida']))
# Agregamos esta nueva información como columna adicional de nuestro shape original de radios.
radios['partidas'] = radios['RADIO_I'].map(d)
# Hacemos una primera visualización de este resultado.
fig = plt.figure(figsize=(15,7))
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)
radios.to_crs(wkt).plot(ax=ax1, column='partidas', scheme='Natural_Breaks', k=3,
linewidth=0.2, edgecolor='black', cmap='Reds_r', legend=True)
radios.to_crs(wkt).plot(ax=ax2, column='partidas', scheme='Quantiles', k=3,
linewidth=0.2, edgecolor='black', cmap='Blues_r', legend=True)
ax1.set_axis_off()
ax2.set_axis_off()
Vamos a introducir otro tema que retomaremos en el notebook siguiente: esquemas de clasificación. Acá simplemente vamos a comparar dos distintos para detectar si existe algún tipo de diferencia. Como no se pueden apreciar contrastes importantes sólo nos detendremos brevemente en la última distribución. Los quantiles.
# reemplazamos los casos sin valor por 0 (en nuestro join espacial muchos polígonos no coincidían con ningún punto)
radios['partidas'].fillna(0,inplace=True)
# calculamos nuestros bins de corte (quantiles).
qcut = list(radios['partidas'].quantile([0, 0.25, 0.5, 0.75, 1]))
# Esto nos devuelve cinco intervalos, o quintiles que dividen a nuestra población en cinco partes iguales
pd.qcut(radios['partidas'], 5, labels=False).value_counts()
Qué significa esto? Básicamente que el primer 20% de nuestros radios censales no tendrá fachadas o partidas dentro del área delimitada por los bordes del polígono. El segundo 20%, es decir entre el 20 y el 40% de los casos tendrá un máximo de 2 partidas. Del 40 al 60 será entre 2 y 5, del 60 al 80% entre 5 y 10 y el último quintil (o 20% de nuestros casos) trendrá entre 10 y 75 partidas. Es decir, el más disperso. Pero como dijimos, este tema lo retomaremos en nuestra próxima clase cuando revisemos algunas clases y métodos de la librería mapclassify.
# nuestros bins de corte
qcut
Ahora, vamos a utilizar la clase Choropleth de Folium para plotear nuestro agrupamiento por radio. Presten atención a la lógica general. Esta no difiere de los otros mapas que hemos visto. Primero agregar el mapa con .Map()(para lo que dejamos seteado un camino que siempre nos devuelva el centro de nuestro gdf a partir de los ceontroides de las latitudes y longitudes). Y luego, ir creando los distintos objetos que iremos agregando como capas superpuestas.
# creamos el mapa
gdf = radios.copy()
centroide = gdf.geometry.centroid
coordenadas = [centroide.y.mean(), centroide.x.mean()]
layer = folium.Map(location=coordenadas,
zoom_start=11,
height = 500,
control_scale=True)
# esto nos permitiría seleccionar entre distintos tiles de base
tiles = ['stamenwatercolor', 'cartodbpositron', 'openstreetmap', 'stamenterrain']
for tile in tiles:
folium.TileLayer(tile).add_to(layer)
# instanciamos la coropleta
fachadas = folium.Choropleth(name='Radios censales',
geo_data=gdf,
data=gdf[[c for c in gdf.columns if c!='geometry']],
columns=['RADIO_I','partidas'],
key_on='properties.RADIO_I',
fill_color='Reds_r',
bins = qcut,
fill_opacity=0.6,
line_opacity=0.2,
highlight=True,
legend_name='Cantidad de partidas por radio censal (CNPHV-2010)',
smooth_factor=1).add_to(layer)
# creamos un tooltip para hacer un hover en nuestros polígonos
pol_hover = folium.features.GeoJsonTooltip(fields=['BARRIO','RADIO_I','partidas'],
aliases=['Barrio:','Radio_id:','Partidas:'])
# Y lo agregamos al layer inicial
folium.GeoJson(gdf,
style_function=lambda x: {"weight":0.35,
'color':'white',
'fillOpacity':0.0},
smooth_factor=2.0,
tooltip=pol_hover).add_to(fachadas);
layer